function [phi J xx] = keffmg_meshcenter
disp('---------------------------------------------')
disp('--------- 1D mesh center diffusion ----------')
disp('---------------------------------------------')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  multi-group diffusion solver for slab geometry,  j.roberts            %
%  this solver uses mesh centered differencing                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BOUNDARY CODITIONS                                                     %
%   0 == vacuum condition (zero incoming current using partial current)  %
%   1 == reflective condition                                            %
%   2 == incident source (if source is zero, same as vacuum)             %
% SOURCE OPTIONS                                                         %
%   0 == fixed-source (volume source)                                    %
%   1 == eigenvalue  (solves via power iteration)                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%clear; %clc;

bcL = 1;        % left boundary condition (do not mix 0 with 1 or 2!)
bcR = 0;        % right boundary condition
it = 1;         % source option

% material and cross-section parameters
% cross-section format:
%    for all materials
%        for all groups
%           read difc(m,g) siga(m,g) nusi(m,g) sigsg1->gi sigsg2->gi ...
% we limit the data only to downscattering

% Simple 1-D problem
numg = 2; % number of groups
numm = 9; % number of materials
file = 'xsec_2g';
[difc,siga,nusi,sigs] = rdxsec(file,numm,numg);

% mesh information
coarse  =   (0:2:20)*10;
fine    =   20*ones(length(coarse)-1,1);
% material placement
matl    =   [  9   2 1 8 4 7 5 3 6  9 ];
% source in each coarse mesh
source    = [ 1 ];       
% boundar source  [left g1 g2 ... ; right g1 ...]
bsrc     = [ 1 
             0 ];   
          
numc   = length(coarse)-1;   % number of course meshes
n      = sum(fine);          % total number of fine meshes

% initialize the finemesh data values
dd = zeros(n,numg);      % diffusion coefficient
aa = zeros(n,numg);      % absorption cross-section
nf = zeros(n,numg);      % nu*fission cross-section
sc = zeros(n,numg,numg); % 
ss = zeros(n,numg);
dx = zeros(n,1);

j = 0;
for i = 1:numc
    dx( (j+1):(j+fine(i)) ) = (coarse(i+1)-coarse(i))/fine(i);
    for g = 1:numg
	    dd( (j+1):(j+fine(i)), g ) = difc( matl(i), g  ); 
	    aa( (j+1):(j+fine(i)), g ) = siga( matl(i), g  );    
	    nf( (j+1):(j+fine(i)), g ) = dx((j+1):(j+fine(i)))*nusi( matl(i), g  ); 
        for gg = 1:numg
            sc( (j+1):(j+fine(i)), g, gg ) = dx((j+1):(j+fine(i)))*sigs( matl(i), g, gg  );
        end
        if it == 0
        %    ss( (j+1):(j+fine(i)), g  ) = dx((j+1):(j+fine(i)))*source(i, g ); 
        end
    end
    j = sum(fine(1:i));
end




if sum( dx >= sqrt(dd(:,end)./aa(:,end)) ) > 0
    disp('meshes are too large!') % want dx < L = sqrt(DiffCoef/SigA)
end

%%%%%%%%%%%%%%%%%%%% COEFFICIENT MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% build the coefficient matrix (just diagonals)
AD = zeros(n,numg);
AL = zeros(n-1,numg);
AU = zeros(n-1,numg);
for g = 1:numg
    % left boundary
    if bcL==1
        b = 1;
    else
        b = 0;
    end
    AU(1,g)   = -2*dd(1,g)*dd(2,g)/(dd(1,g)*dx(2)+dd(2,g)*dx(1));
    AD(1,g)   = -AU(1,g) + 2*dd(1,g)*(1-b)/...
                    (4*dd(1,g)*(1+b)+dx(1)*(1-b)) + dx(1)*aa(1,g); 
    % left boundary source
    if bcL == 2
        ss(1,g) = ss(1,g) + 8*dd(1,g)*bsrc(1,g)/(4*dd(1,g)+dx(1));
    end  

 
    % right boundary
    if bcR==1
        b = 1;
    else
        b = 0;
    end
    AL(n-1,g) = -2*dd(n,g)*dd(n-1,g)/(dd(n,g)*dx(n-1)+dd(n-1,g)*dx(n));
    AD(n,g)   = -AL(n-1,g) + 2*dd(n,g)*(1-b)/...
                    (4*dd(n,g)*(1+b)+dx(n)*(1-b)) + dx(n)*aa(n,g);               
    % right boundary source
    if bcR == 2
        ss(end,g) = ss(end,g) + 4*dd(end,g)*bsrc(end,g)/(4*dd(end,g)+dx(end));  
    end          
                
    for i = 2:n-1  % interior
        AL(i-1,g) = -2*dd(i-1,g)*dd(i,g)/(dd(i-1,g)*dx(i)+dd(i,g)*dx(i-1));
        AU(i,g)   = -2*dd(i+1,g)*dd(i,g)/(dd(i+1,g)*dx(i)+dd(i,g)*dx(i+1));   
        AD(i,g)   = -( AL(i-1,g) + AU(i,g) ) + dx(i,1)*aa(i,g);
    end
end

if it == 1
    nusig = zeros(n,numg);
else
    s = zeros(n,numg);
end
sct = zeros(n,numg,numg);

if it == 1
    sindx = find( sum(nusig(:,:)') ); % if no nusig, so source  
    s = 0.4*ones(n,1).*( sum(nusig(:,:)' > 0 ) )';  
    s = 0.01 * s / max(s);
end

%------------- make the x-vector for plotting------------------------
x = zeros(length(dx)+1,1);
for i = 2:1:length(dx)+1
	x(i) = x(i-1,1) + dx(i-1,1);
end

phi     = zeros(n,g);
scsrc   = zeros(n,1);
if it == 1
    s = dx.*ones(n,1).*(nf(:,2)); 
end
s = s/norm(s);
if it == 1  % eigenvalue search -------------------------------------
    [k s phi] = keffit( AU, AL, AD, dx, nf, sc, s, 1e-5, 1e-5, 100);
else        % fixed source, probably want function ------------------
    for g = 1:numg
        % compute scattering source
        for k = 1:n
            for gg = 1:g-1 % only down scatter
                scsrc(k,1) = scsrc(k,1) + sc(k,g,gg)*phi(k,gg);
            end
        end
        % solve for phi_g
        phi(:,g) = tridiag(AU(:,g),AL(:,g),AD(:,g),ss(:,g)+scsrc(:,1));
    end  
end

% J(x) ~ - D(x)*(d/dx)phi(x)
% phi_prime_edge = zeros(length(phi,1)+1,1);
% phi_prime_edge(1)    = 2/(4*dd(1) + dx(1)) * phi(1);
% phi_prime_edge(end)  = 2/(-4*dd(end) + dx(end)) * phi(end);
% phi_prime_edge(2:end-1) =  2*(phi(2:end)-phi(1:end-1)) .* dd(2:end) ./ ( dx(1)*dd(1:end-1)+dx(1)*dd(2:end));
% phi_prime = 0.5*(phi_prime_edge(2:end)+phi_prime_edge(1:end-1));
J = 0;%*(phi + 0.5*-dd.*phi_prime);

xx = 0.5*(x(1:end-1)+x(2:end));

% figure(1)
% if numg==1
%     plot(xx,phi(:,1),'k',xx,J,'b--','LineWidth',2)
%     axis([min(x) max(x) 1.25*min(J) 1.25*max(max(phi))]);
%     xlabel('{x(cm)}')
%     ylabel('\phi(x) or J(x) {n/cm^2-s}')
%     legend('\phi(x)','J(x)')
%     grid on
% elseif numg==1
%     plot(xx,phi(:,1),'k-.','LineWidth',2)
%     axis([min(x) max(x) 1.25*min(J) 1.25*max(max(phi))]);
%     xlabel('{x(cm)}')
%     ylabel('\phi(x)')
% end


%max(phi(:,1)./phi(:,2))

end

function [k,s,phi,iter] = keffit(c,b,a,dxx,nusig,sct,s,epsK,epsS,itmx )

k       = 1.0;    % guess the initial keff
errK    = 1; 
errS    = 1; 
iter    = 0; 
sindx   = find(s);
n       = length(dxx);
numg    = length(nusig(1,:));
phi     = zeros(n,numg); 

while   ( ( errK>epsK ) || (errS>epsS) ) && iter < itmx
    % call tridiagonal solver --------------------------
    %s = s/sum(s);
    phi(:,1)    = tridiag(b(:,1),c(:,1),a(:,1),s/k);
    if numg > 1                     % successive groups
        for g = 2:numg
            scsrc = zeros(n,1);
            for gg = 1:g-1        
            	 scsrc(:,1) = scsrc(:,1) + sct(:,g,gg).*phi(:,gg);
            end 
            % solve for phi_g
            phi(:,g)  = tridiag(c(:,g),b(:,g),a(:,g),scsrc(:,1));
        end
    end
    sold = s; 
    kold = k;   
    s(:) = 0;       % reset fission source
    for g = 1:numg  % compute fission source
      s(:) = s(:) + nusig(:,g).*phi(:,g);
    end  
    % compute keff
    k = sum(s.*dxx)*kold/sum(sold(sindx).*dxx(sindx)); 
    errS = max( abs((s(sindx)-sold(sindx))./s(sindx)) );
    errK = abs( (k-kold)/k );
    iter = iter + 1;  % number of iterations
end
k
disp(['*** final keff estimate: ', num2str(k),' ***'])
disp(['residual errors: errS = ',num2str(errS),' and errK = ',num2str(errK)])
disp(['in ',num2str(iter),' iterations'])

end

